Data Diary: Scenarios

Background

To develop the scenarios for nitrogen loss in Rock Your Watershed, we used water monitoring data from 4 small catchments in Iowa. In addition, we used climate data from each catchment. Both data sets are at daily intervals from 2015 through 2019.

Import

The data from all four catchments (nitrogen concentration, flow, and precip) was compiled into one XLSX file, with catchment’s individual data on separate sheets. To import it, we made a function that goes one tab at a time and renames some of the columns for easier handling. In addition, it calculates the N loading (c*Q) and creates a seasonal value based on the month for future aggregations.

site_fc <- function(number){
  
  # one catchment per sheet
  c1 <- read_excel("Data/Peiyu/KS_RS_LP_WW_for_Ignacio.xlsx", sheet = number)
  
  # rename columns and remove empty ones
  c1 <- select(c1, 
               date = Date,
               flow_m3_day = `Inflow (m3/day)`,
               c_mg_L = `[NO3] in (mg/L)`,
               precip_mm = `Precipitation(mm)`)
  
  # now create columns for:
  # N load, breakdown of date/season for future aggregations
  
  c1 <- c1 %>% mutate(load_g_day = round(c1$flow_m3_day * c1$c_mg_L), 
                      site = number,
                      year = year(date),
                      month = month(date),
                      day = day(date),
                      season = case_when(
                        month %in% c(12, 1, 2) ~ "1_Winter",
                        month %in% c(3:5) ~ "2_Spring",
                        month %in% c(6:8) ~ "3_Summer",
                        month %in% c(9:11) ~ "4_Fall"
                      ))
  
  return(c1)
  
}

This function was iterated 4 times, once for each site. Catchments were also assigned a site number (1-4) and then stored on a list. This can be used as a small “database” to access each site’s dataframe by index number.

site_list <- lapply(1:4, site_fc)

Extreme Precipitation

For this analysis, we defined extreme precipitation as the top 10% of precipitation days by rain amount. To avoid skewed data, the percentile threshold is defined by comparing precip values on a monthly basis from all years of the same site.

The following function makes this breakdown and assigns binary values to each day:

• 0: nonextreme
• 1: extreme

The function is applied to each site individually and then all the data is stored in a new list.

Note: days with precip below 1 mm are discarded for the threshold calculation.

## the function inputs are the site number and percentile threshold

extreme_precip_fc <- function(number, perc){
  
  m1 <- as.data.frame(site_list[number]) %>%
    filter(month == 1) 
  
  m1 <- m1 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m1, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m1, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m2 <- as.data.frame(site_list[number]) %>%
    filter(month == 2) 
  
  m2 <- m2 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m2, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m2, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m3 <- as.data.frame(site_list[number]) %>%
    filter(month == 3) 
  
  m3 <- m3 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m3, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m3, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m4 <- as.data.frame(site_list[number]) %>%
    filter(month == 4) 
  
  m4 <- m4 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m4, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m4, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m5 <- as.data.frame(site_list[number]) %>%
    filter(month == 5) 
  
  m5 <- m5 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m5, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m5, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m6 <- as.data.frame(site_list[number]) %>%
    filter(month == 6) 
  
  m6 <- m6 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m6, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m6, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m7 <- as.data.frame(site_list[number]) %>%
    filter(month == 7) 
  
  m7 <- m7 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m7, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m7, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m8 <- as.data.frame(site_list[number]) %>%
    filter(month == 8) 
  
  m8 <- m8 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m8, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m8, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m9 <- as.data.frame(site_list[number]) %>%
    filter(month == 9) 
  
  m9 <- m9 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m9, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m9, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m10 <- as.data.frame(site_list[number]) %>%
    filter(month == 10) 
  
  m10 <- m10 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m10, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m10, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m11 <- as.data.frame(site_list[number]) %>%
    filter(month == 11) 
  
  m11 <- m11 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m11, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m11, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  m12 <- as.data.frame(site_list[number]) %>%
    filter(month == 12) 
  
  m12 <- m12 %>% 
    mutate(extreme_day = case_when(
      precip_mm >= quantile(with(m12, precip_mm[precip_mm >= 1]), perc) ~ 1,
      precip_mm < quantile(with(m12, precip_mm[precip_mm >= 1]), perc) ~ 0
    ))
  
  
  m <- rbind(m1, m2, m3, m4, m5, m6,
             m7, m8, m9, m10, m11, m12)
  
  rm(m1, m2, m3, m4, m5, m6,
     m7, m8, m9, m10, m11, m12)
  
  return(m)
  
}

## apply the function to all sites and compile the results
extreme_list <- lapply(1:4, extreme_precip_fc, 0.9) 

Scenario Levels

Once the extreme precipitation values are assigned, data was aggregated by catchment-season. On this aggregation, we calculated the sum of precipitation (mm), loading (g) and number of extreme precipitation days.

For the first layer of climate scenarios, we defined the precipitation levels like this:

• High precip: 20% higher than the mean value for that same season type
• Normal precip: +-20% of the mean value for that same season type
• Low precip: 20% lower than the mean value for that same season type

The threshold values are once again calculated by comparing data points from the same season to avoid bias (i.e. spring/summer values tend to be higher than winter). Initially, we used 10% above mean values but that shifted the vast majority of data points to either high or low. Using 20%, there is a better distribution of data points across the three levels.

For the second level of climate scenarios, we defined the extreme precipitation level like this:

• More extreme: a season with a number of extreme days greater than the median value for that precipitation level (High/Normal/Low)
• Less extreme: a season with a number of extreme days less than or equal to the median value for that precipitation level (High/Normal/Low)

The following function performs these calculations for each site individually and then the results are compiled.

precip_season_fc <- function(number){
  
  season1 <- as.data.frame(extreme_list[number]) %>%
    group_by(year, season) %>%
    summarise(n_days = sum(extreme_day),
              precip_mm = sum(precip_mm),
              load_g = sum(load_g_day),
              flow_m3 = sum(flow_m3_day)) %>%
    mutate(site = number)
  
  # precip level determined like extreme day, comparing all same season type of same site
  
  s1 <- season1 %>%
    filter(season == "1_Winter") 
  
  s1 <- s1 %>% 
    mutate(precip_level = case_when(
      precip_mm >= mean(s1$precip_mm)*1.2 ~ "3_High",
      precip_mm <= mean(s1$precip_mm)*0.8 ~ "1_Low",
      TRUE ~ "2_Normal"
    ))
  
  s2 <- season1 %>%
    filter(season == "2_Spring") 
  
  s2 <- s2 %>% 
    mutate(precip_level = case_when(
      precip_mm >= mean(s2$precip_mm)*1.2 ~ "3_High",
      precip_mm <= mean(s2$precip_mm)*0.8 ~ "1_Low",
      TRUE ~ "2_Normal"
    ))
  
  s3 <- season1 %>%
    filter(season == "3_Summer") 
  
  s3 <- s3 %>% 
    mutate(precip_level = case_when(
      precip_mm >= mean(s3$precip_mm)*1.2 ~ "3_High",
      precip_mm <= mean(s3$precip_mm)*0.8 ~ "1_Low",
      TRUE ~ "2_Normal"
    ))
  
  s4 <- season1 %>%
    filter(season == "4_Fall") 
  
  s4 <- s4 %>% 
    mutate(precip_level = case_when(
      precip_mm >= mean(s4$precip_mm)*1.2 ~ "3_High",
      precip_mm <= mean(s4$precip_mm)*0.8 ~ "1_Low",
      TRUE ~ "2_Normal"
    ))
  
  season2 <- rbind(s1, s2, s3, s4)
  
  # separate precip_level by precip levels (low, normal, high)
  # then assign binary values for high and low # of extreme days
  
  season_1_Low <- season2 %>% filter(precip_level == "1_Low")
  season_2_Normal <- season2 %>% filter(precip_level == "2_Normal")
  season_3_High <- season2 %>% filter(precip_level == "3_High")
  
  season_1_Low <- season_1_Low %>% mutate(extreme_level = case_when(
    n_days <= mean(season_1_Low$n_days) ~ "1_Less_Extreme",
    n_days > mean(season_1_Low$n_days) ~ "2_More_Extreme"
  ))
  
  season_2_Normal <- season_2_Normal %>% mutate(extreme_level = case_when(
    n_days <= mean(season_2_Normal$n_days) ~ "1_Less_Extreme",
    n_days > mean(season_2_Normal$n_days) ~ "2_More_Extreme"
  ))
  
  season_3_High <- season_3_High %>% mutate(extreme_level = case_when(
    n_days <= mean(season_3_High$n_days) ~ "1_Less_Extreme",
    n_days > mean(season_3_High$n_days) ~ "2_More_Extreme"
  ))
  
  season3 <- rbind(season_1_Low, season_2_Normal, season_3_High)
  
  rm(s1, s2, s3, s4, season_1_Low, season_2_Normal, season_3_High, season2)
  
  return(season3)
  
}


# run an iteration of the function for each site (1-4)

for (i in 1:4) {
  
  name <- paste("precip_season", i, sep = "_")
  
  assign(name, precip_season_fc(i))
  
}

# this time combine the results into one dataframe
precip_season <- rbind(precip_season_1, precip_season_2, precip_season_3, precip_season_4)

Visualization

These are the data points once they are aggregated by season and classified into the first layer of precipitation level:

viz_linear

Here, they are organized by scenario and subscenario. As a reference, the red line shows the mean value for all catchment-season points.

viz_box_load

Results

To implement the 6 nitrogen loss subscenarios into the game, we calculated the percentage change of mean loading for each. The baseline for the percentage change is the mean loading across all catchment-season data points.

table_1_Low <- aggregate(box_1_Low[, 5], list(box_1_Low$extreme_level), FUN=mean) %>%
  mutate(baseline = mean(precip_season$load_g),
         pct_change_load = round((load_g - baseline)/baseline, 2),
         precip_level = "1_Low")

table_2_Normal <- aggregate(box_2_Normal[, 5], list(box_2_Normal$extreme_level), FUN=mean) %>%
  mutate(baseline = mean(precip_season$load_g),
         pct_change_load = round((load_g - baseline)/baseline, 2),
         precip_level = "2_Normal")

table_3_High <- aggregate(box_3_High[, 5], list(box_3_High$extreme_level), FUN=mean) %>%
  mutate(baseline = mean(precip_season$load_g),
         pct_change_load = round((load_g - baseline)/baseline, 2),
         precip_level = "3_High")

results <- rbind(table_1_Low, table_2_Normal, table_3_High)

scenario_name_v <- c("1: Low Precip/Less Extreme", "2: Low Precip/More Extreme", 
                     "3: Normal Precip/Less Extreme", "4: Normal Precip/More Extreme",
                     "5: High Precip/Less Extreme", "6: High Precip/More Extreme")

results <- results %>%
  mutate(scenario_name = scenario_name_v,
         coeff = 1 + pct_change_load) %>%
  select(scenario_name, load_g, baseline, pct_change_load, coeff)

results %>% 
  ggplot(aes(x = scenario_name, y = pct_change_load, fill = pct_change_load)) +
  geom_bar(stat="identity") +
  scale_y_continuous("Load, percentage change") +
  scale_x_discrete("") +
  scale_fill_distiller(palette = "BrBG", limits = c(-1,1)) +
  theme_linedraw()

As you can see in the plot above, both higher the precipitation level and a higher number of extreme precipitation days, can lead to greater nitrogen loss.

We also performed a similar analysis on water flow instead of nitrogen loading. These are our other results: